//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "DGDirichBC.h"

// MOOSE includes
#include "MooseVariableFE.h"

#include "libmesh/numeric_vector.h"
#include "libmesh/utility.h"

// C++ includes
#include <cmath>

registerMooseObject("snakeApp", DGDirichBC);

template <>
InputParameters
validParams<DGDirichBC>()
{
  InputParameters params = validParams<IntegratedBC>();
  params.addRequiredParam<Real>("value", "The value the variable should have on the boundary");
  params.addRequiredParam<Real>("epsilon", "Epsilon");
  params.addRequiredParam<Real>("sigma", "Sigma");
   params.addParam<MaterialPropertyName>(
      "kappa", "kappa", "Property name of the  heat conductivity material property");
  return params;
}

DGDirichBC::DGDirichBC(const InputParameters & parameters)
  : IntegratedBC(parameters),
    _epsilon(getParam<Real>("epsilon")),
    _sigma(getParam<Real>("sigma")),
    _bcval(getParam<Real>("value")),
    _kappa(getMaterialProperty<Real>("kappa"))
{
}

Real
DGDirichBC::computeQpResidual()
{
  const unsigned int elem_b_order = _var.order();
  double h_elem =
      _current_elem->volume() / _current_side_elem->volume() * 1. / Utility::pow<2>(elem_b_order);
  

  Real r = 0;
  r -= (_kappa[_qp] * _grad_u[_qp] * _normals[_qp] * _test[_i][_qp]);
  r += _epsilon * (_u[_qp] - _bcval) * _kappa[_qp] * _grad_test[_i][_qp] * _normals[_qp];
  r += _sigma / h_elem * (_u[_qp] - _bcval) * _kappa[_qp]*_test[_i][_qp];

  return r;
}

Real
DGDirichBC::computeQpJacobian()
{
  const unsigned int elem_b_order = _var.order();
  double h_elem =
      _current_elem->volume() / _current_side_elem->volume() * 1. / Utility::pow<2>(elem_b_order);

  
  Real r = 0;
  r -= (_kappa[_qp] * _grad_phi[_j][_qp] * _normals[_qp] * _test[_i][_qp]);
  r += _epsilon * _phi[_j][_qp] *_kappa[_qp] * _grad_test[_i][_qp] * _normals[_qp];
  r += _sigma / h_elem * _phi[_j][_qp] *_kappa[_qp]* _test[_i][_qp];

  return r;
}
